Finite-Temperature Neel Ordering of Fluctuations in a Plaquette Orbital Model 
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We present a pseudo-spin model which should be experimentally accessible using solid-state devices and, 
being a variation on the compass model, adds to the toolbox for the protection of qubits in the area of quantum 
information. Using Monte Carlo methods, we find for both classical and quantum spins in two and three di- 
mensions Ising type Neel ordering of energy fluctuations at finite temperatures without magnetic order. We also 
readdress the controversy concerning the stability of the ordered state in the presence of quenched impurities 
and present numerical results which are at clear variance with earlier claims in the literature. 
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I. INTRODUCTION 

The prospect of topological quantum computation to imple- 
ment fault tolerant quantum bits has led to considerable inter- 
est in the fieldA over the past years. A particular route within 
this area is the construction of simple (spin) models, as micro- 
scopic models of topological field theories, that allow a direct 
experimental realization. A hallmark is the so-called Kitaev 
model^ which is exactly solvable and known to be imple- 
mentable as well as controllable using optical devices^ Very 
recently, related efforts have been spent to construct similar 
models which can be realized using solid-state techniques. 4 
Guided by a few principles - a degenerate ground-state mani- 
fold and a gap to excited states - the so-called compass model 
(CM), with the Hamiltonian 



(1) 



on an L x L lattice, was proposed as a simple model allow- 
ing for the protection of qubits. 5 The CM is realizable by a 
particular arrangement of Josephson-junction devices and can 
be described by a Z2 Chern-Simons topological quantum field 
theory; 5 Extensions of (Q]) to global interactions possess even 
better fault tolerant properties^ 

Originally introduced as an orbital model for Mott 
insulators^ research into the actual CM has been pushed 
by several groups, which established degenerate ground-state 
properties^ first-order quantum phase transitions^ relation to 



— and the existence of directional 
argues that quantum spins support 



p + ip superconductivit 
orderi 11 ' 12,13 Reference [T 
a resistivity of the ordered phase towards quenched disor- 
der, which is in sharp contrast to classical degrees of free- 
dom for which the ordered phase vanishes rapidly with in- 
creasing disorder. By a detailed Monte Carlo (MC) studyi^ 
of the quantum and classical CM, we have recently shown 
that this conclusion, however, needs further support as the 
CM shows unusual and extremely slowly converging finite- 
size scaling (FSS) properties on periodic lattices, which were 
used in Ref. 12;. Recently, an interesting extension of (fl~|i to 
include a magnetic field term hS was performed, leading to 
thermal canting of spin orderi 

In the search for other fundamental spin models and to gain 
further insights into the field around the Kitaev model, we pro- 



pose here a different - geometric - modification of the CM 
and concisely report on its intriguing physics. Our main result 
is the establishment of an interesting ordering that can be de- 
scribed by a crystallization and modulation of local energy 
contributions but which lacks conventional magnetic order. 
We show that the proposed model falls into the Ising univer- 
sality class and that it possesses well behaved FSS properties 
in contrast to the CM. In the last part of this paper, we use this 
advantage to investigate the influence of (weak) quenched dis- 
order in form of random vacancies to study their influence on 
the nature of the phase transition. We show that long-range 
order is completely lost already for very weak impurity con- 
centrations. 



H. THE MODEL 

The plaquette orbital model (POM) is defined by the Hamil- 
tonian 



'HpOM — J A ^2 ^iSj + Jb 
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<»>J>S 



QZ H2 



(2) 



where S x and S z are components of a two-component spin 
S, which can represent both classical and quantum degrees 
of freedom. In the latter case S x and S z are represented by 
the usual 5 = 1/2 Pauli matrices while in the classical case 
they denote projections of a continuous spin parameterized by 
an angle 9 on the unit sphere. The bonds (i,j)A, {i,J)B on 
sub-lattices A and B are arranged as depicted in Fig. [TJ The 
coupling strengths J a and Jb are in principle arbitrary. Here, 
we are interested in the isotropic case J a = Jb = J = 1. The 
sign of J has no relevance since it can be transformed away 
on bipartite lattices^ With N = L d we denote the number of 
spins on a cubic lattice of linear extension L and dimension 
d. It should be emphasized that in contrast to the CM in three 
dimensions (3D) or the Kitaev model in two dimensions (2D), 
quantum MC investigations of the POM can easily be done 
also in 3D since there is no sign problem. 

Note that the Hamiltonian (O is Z2 symmetric under ex- 
change of sub-lattices A and B and spin indices x and z. De- 
fine further four-spin operators 



P, 



C2 02 QZ QZ 



(3) 



2 







B 








A 












B 

sis) 








SfS] 
1 j 

A 


( 1 




( 













FIG. 1: (color online). Illustration of the POM lattice. The blue 
(thick) bonds indicate S x S x terms while the red (dashed) bonds 
stand for S Z S Z links. The lattice is closely related to a checker- 
board. The generalization to the 3D POM is straightforward with 
sub-lattices A and B being cubes rather than plaquettes. 
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and 

Qi = Sf Sf +ex Sf +ey Sf +ey +Cx , (4) 

with and e y being unit vectors on the lattice and r (/) 
a site pointing to the lower left corner of A (B) plaque- 
ttes. We can show that [HpoM,P r } = [Wpom, Qi] = 
are local symmetries. Hence, on A plaquettes the operations 
Sf — > —Sf , Sf — > Sf are symmetries with the analog ex- 
pression on B plaquettes reading Sf — > S*f , Sf — > — SJ^ For 
any r and L > 2, there further exists at least one index I such 
that [P r , Qi] ^ 0,- This shows that every energy eigenvalue 
of 7Ypom is at least two-fold degenerate. Performing an exact 
diagonalization employing invariant subspaces of the operator 
P (see Ref. |8|), we could obtain all eigenvalues on a N = 4 x 4 
cluster confirming this conclusion. The POM hence possesses 
the same behavior as the CM in this regard. Whether the ex- 
citation gap persists in the infinite-volume limit remains to be 
investigated. The 3D extension of the model is obvious. Ev- 
ery plaquette becomes a cube, otherwise all arguments stay 
the same. 

We now turn to a stochastic investigation of the model us- 
ing established MC methods ranging from Metropolis sam- 
pling for classical variables to the quantum stochastic series 
expansion (SSE)4&ii Key to succesfully simulate model (O 
on large lattices is the use of parallel tempering methods to 
avoid barriers and reduce autocorrelation times. Details of 
our approach can be found in Ref. [TH 



III. RESULTS 

A. Neel ordering 

The exchange symmetry (A <^> B) provides the possibility 
for spontaneous symmetry breaking. To see if this exchange 
symmetry is broken for some temperature T < T c a suitable 



FIG. 2: (color online). Snapshots of the plaquette energy distribution 
on a N = 40 x 40 lattice taken at T = 1.0 and T = 0.02, respec- 
tively. In the high-temperature regime (a) the system is disordered 
and a symmetry breaking has occurred in (b) for T less than a crit- 
ical temperature T c . Darker regions (color) represent lower energy. 
Figures (c) and (d) are snapshots of configurations in spin space cor- 
responding to the high- and low-temperature phases for L = 10. No 
evident magnetic order is seen. The circles in (c) and (d) signify local 
energy density where darker means larger negative energy. Spins are 
color (gray) coded to make their direction more apparent. 

order parameter should be defined. Let us just consider energy 
fluctuation on the two sub-lattices. Then, following Refs. [l2l 
[l3lfT8l we define the quantity 

D=^\E A -E B \, (5) 

where Ea = J a j\ A SfSJ is the energy contribution on 
sub-lattice A with the obvious relation for Eb- The suitabil- 
ity of that quantity can be seen in Fig. |2|a,b) which shows 
snapshots of the energy distribution at high and low tempera- 
tures. We can clearly observe a phase transition and the low- 
temperature phase can be described as a crystalline Neel state 
expressing up-down energy modulations. Further, this state is 
entirely described by energy fluctuations as there is no sign of 
long-ranged magnetic order seen in Fig.[2|c,d), as expectediJA 
A quantity directly probing a crystalline state as in Fig. |2b) 
is for example a plaquette structure factor defined by 

jv 

S p i = (l/A0£(-l) a ' + *Sp(0)£p(r), (6) 

r 

where 

Sp(r) = SVSV+e^ + S r J rey S r ^ e _ c + ey + S r S r -\-e v 

~\~ &T+e x Sr+e<t+ey (7) 
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is a four-site spin operator on a plaquette. The sign of 
(— l) x ' ' +Vr alternates between A and B plaquettes. We will 
show below that S v \ is indeed an order parameter. 



B. Monte Carlo simulations 

Having visualized the onset of crystalline order, we now 
turn to a more careful discussion of the phase transition from 
comprehensive MC runs. Finite lattices are studied using two 
types of boundary conditions. In principle, we do not ex- 
pect unwanted excitation as in the CM which spoil the FSS 
on periodic lattices^ but we also study open boundary condi- 
tions to gain further confidence in our results. Open boundary 
conditions have the additional advantage that they prefer one 
sub-lattice over the other thus possibly stabilizing the ordered 
phase from the surface. In the latter case the | . . . | in Eq. (0 
can be also left away from the definition of the order param- 
eter. In 2D, simulations were performed of both the classical 
and quantum cases for various lattice sizes L = 10, ... , 96, 
which proved to be sufficiently large. Our analysis to obtain 
the critical temperatures is based on D, rather than S p \ be- 
cause it is easier measured and is less susceptible to statistical 
noise. Detection of the phase boundary proceeds, as usual, by 
considering the susceptibility 



N((D- 2 )-(Df), 
or the Binder parameter defined as 

(D 4 )/ (3(DY) 



\ 



B = 1 



(8) 



(9) 



Figure [3{a) summarizes data obtained for the classical model 
for the order parameter D as well as the structure factor 5 p i 
(in inset). Both quantities clearly numerically establish exis- 
tence of crystalline order. In Fig.[3b), data for B and the sus- 
ceptibility x (in inset) is given which suggests a second-order 
phase transition at T c;c \ = 0.0855(4) from the crossings of the 
Binder parameter. The value of B(T C ) = 0.610(5) at the crit- 
ical temperature is consistent with the usual 2D Ising value on 
the torus topologyJ2ia£ We further investigate the critical tem- 
perature by a FSS analysis of the maxima of the susceptibility 
X- By fitting the corresponding data in Fig. [^c) to the usual 
ansatz 



r roa x(i) =T C + aL 



bL~ 



(10) 



we obtain T c;cl = 0.0860(2) and T c;qu = 0.048(1) for the 
classical and quantum cases, respectively. Here, we assume 
the 2D Ising exponent v = J^Dismg = 1 (justified by the 
fit quality and independent fits to the slope of the Binder pa- 
rameters), and the effective correction term ~ L~ w is used 
only in fits for open boundary conditions. Since these crit- 
ical values are consistently obtained for open and periodic 
boundary conditions, we arrive at the important conclusion 
that FSS in the POM is well behaved, in clear distinction to 
the CM. It is then instructive to compare those values with the 
related critical temperatures of the directional-ordering tran- 
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FIG. 3: (color online). Results from classical and quantum MC simu- 
lations for the 2D POM. (a) The order parameter D for various lattice 
sizes L = 10 to L — 96 for the classical model employing periodic 
boundary conditions. The inset shows the plaquette structure factor 
S p i for L = 20, L = 32, and L = 42. The order parameter indi- 
cates a clear crystallization effect for T < 0.086. (b) The classical 
Binder parameter B close to the transition temperature supporting 
a second-order phase transition. The inset shows the susceptibility 
X on a logarithmic scale for L > 16. Generally, steeper curves in 
(a) and (b) correspond to larger lattice sizes, (c) The critical tem- 
peratures T c from FSS of the maxima locations of the susceptibility 
for both the classical and quantum cases and two different boundary 
conditions (open = obc, periodic = pbc). 



and T c:qu = 0.055(1), which leads us to the conclusion that 
the geometric variation from Hamiltonian ([T]i to (O results in 
a drastic reduction of T c by 42% for classical fields vs only 
13% for the quantum case. 

In order to obtain further clarity on the type of the tran- 
sition, we go one step further and study the model in three 
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FIG. 4: (color online). FSS of the order parameter D for a 2D sys- 
tem and for the slope sb of the Binder parameter for a 3D system. 
Both quantities are calculated at the critical temperatures and yield 
the expected Ising exponents. 

dimensions for various lattice sizes L = 8, . . . , 32. Without 
showing details but by performing the same simulations and 
FSS analysis as before, we obtain a clear signal for a long- 
ranged cube-ordered state below the transition temperatures 
T c 3 ° = 0.365(1) and T c 3 £ u = 0.180(5) for classical and 
quantum degrees of freedom, respectively. Interestingly, the 
increase in T c (compared to the 2D transition temperatures) is 
larger for the classical model. 

The proposition that the transition is of Ising type (as sug- 
gested by the symmetry and Binder parameter in Fig. 13b)) 
should be reflected in the critical exponents. In Fig. |4]we se- 
lect two quantities to address this question, namely the order 
parameter D for the 2D case and the slope of the Binder 
parameter at the critical point for the 3D case. By performing 
fits to the ansatz D ~ L~^l v and sb ~ L 1 ^, we obtain in 
the classical cases (/?/V) 2 d = 0.124(2) and v w = 0.62(2), 
which are in excellent agreement with the theoretical value 
(/3/^)2D = 1/8 and high-precision literature^ on z/3D- Fig- 
ure|4]shows this scaling versus the lattice size of both observ- 
ables at the critical temperatures in a log-log plot. The good 
quality of our data is apparent. This establishes that the tran- 
sition in the POM has Ising exponents for both 2D and 3D, 
and - more importantly - that the energetic quantity D really 
scales like a (magnetic) order parameter. For the quantum 
case, this analysis is not so easy but our data is consistent with 
this conclusion. 



C. Dilution effects in the POM 

As the preceding analysis indicates that we have good FSS 
behavior in the present model, it is well suited to re-address 
the important question of impurity effects^ in orbital mod- 
els. To this end, we employ periodic boundary conditions 
and define a fraction x of quenched impurity sites, i.e., we 
remove each spin with probability x. Following Ref. [l~2l our 
objective is to study the quantity g(x) = T c (x) /T c (0) to learn 
about the degree of stability of the Neel-ordered phase against 
dilution disturbances. With T c (x) one refers to the critical 
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FIG. 5: Plot of (a) the disorder averaged (order) parameter [_D] av and 
(b) the disorder averaged Binder ratio [B] av versus temperature T 
for x = 1% site dilution in the classical case. The Binder parameter 
does not show any signature of a remaining phase transition for sizes 
up to L = 128 and temperatures down to T = 0.03. The critical 
temperature in the clean case is indicated by the arrow. 



temperature of the phase transition with dilution x. In or- 
der to access T c (x), we performed simulations on lattice sizes 
L = 20,32,42,64, 128 (classical case) and L = 20, 32 (quan- 
tum case), where we generated and studied 100 — 200 differ- 
ent disorder realizations, respectively. The quantities [-D] av , 
[x]av, and [B] av denote the disorder averaged values of the 
respective quantities defined above. 

Studying just the (peaks of the) averaged susceptibility 
or the order parameter on moderate lattices sizes (L = 
20,32,42), one could deduce values for g(x) which are of 
the order of g(x) w 0.70 for the classical case vs g(x) w 0.82 
for the quantum model in case of weak impurity concentra- 
tion with x = 0.01. These values seem to be in qualitative 
agreement with earlier claims and in support of the conclu- 
sion of Ref. [12 that quantum fluctuations make the ordered 
phase more robust. 

However, our simulations on bigger lattices for x = 0.01 
reveal that the transition is, in fact, vanishing in the thermody- 
namic limit, implying g(x) = 0, V.t > 0. This conclusion can 
be drawn for example from the data shown in Fig. [5] While 
the order-parameter [D] av seems to indicate some crossover 
which gets weaker for larger L, the Binder parameter [B] av 
clearly shows no crossing, i.e., no sign of critical behavior. 
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Remarkably, quite large lattice sizes L « 128 are needed to 
see this. These data obviously rule out the presence of a true 
phase transition in the vacancy diluted classical POM. 

This situation is further reminiscent of the ordinary 2D 
Ising model subjected to a random field (at a fraction x of the 
sites) which is known not to exhibit a finite-temperature phase 
transition from theoretical 2 * and numerical worksS. While 
in the Ising model the random field destroys the Z2 up-down 
spin symmetry, dilution in the POM breaks the A — B pla- 
quette symmetry. Such random local symmetry breaking is 
also known to destroy long-range order in, e.g., 2D antiferro- 
magnetic Ising models with nearest and next-nearest neighbor 

OA 

interactions^ 

By the same argument it is clear that there is no phase tran- 
sition in the 2D classical CM at any finite dilution x and the 
statements and conclusions of Ref.[l2|are therefore at variance 
with the findings in this work. We also see no argument why 
the quantum CM should behave differently in this respect and 
suspect that quantum effects merely increase the stability of 
the low-temperature state on finite and small clusters of spins 
— an observation that might still be useful for applications. 

IV. CONCLUSION 

We have introduced and investigated a plaquette orbital 
model related to the Kitaev model and the CM. The present 



work establishes that this model exhibits an unconventional 
finite-temperature phase transition from a disordered to a 
Neel-ordered state in the energy distribution. It thus displays 
antiferromagnetic features without possessing magnetic or- 
der. By symmetry arguments and MC simulations, the crit- 
ical exponents were determined to be in the Ising universality 
class. The geometric variation from the CM to the plaquette 
orbital model results in a substantial lowering of the order- 
ing temperature for the classical model, which is not the case 
for quantum spins. Our subsequent analysis of the POM in 
the presence of impurities shows that long-range ordering is 
lost for (any) weak disorder concentration. This finding also 
sheds concluding light on the somewhat controversial issue 
concerning the effect of impurities on the ordered phase in 
the compass modelji 2 *^ A more detailed analysis of the na- 
ture of the ground states for arbitrary J a, Jb and its quan- 
tum phase transitions 9,25 would be an interesting continuation 
of this work as would be a thorough investigation of the ex- 
citations in the POM. Apart from its possible relevance for 
protected qubits, the present model should also be of imme- 
diate interest for the physics of orbital models in relation to 
transition-metal oxides^ 
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